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Abstract 



The purpose of these lectures is to give a short introduction into a very vast field of 
numerical simulations for cosmological applications. I focus on major features of the simu- 



lations: the equations, main numerical techniques, effects of resolution, and methods of halo 
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identification. 
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OO . 1- Introduction 
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Numerical simulations play a very significant role in cosmology. It all started in 70s with 
simple N-body problems solved using primitive N-body codes with few hundred particles. 
Later the code which we now call Particle-Particle code (direct summation of all two-body 
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forces) was polished and brought to the state-of-art [1 ] . Already those early efforts brought 
some very valuable fruits. In 1970 Peebles [2] studied collapse of a cloud of particles as a 



model of cluster formation. The model had 300 points initially distributed within a sphere 
with no initial velocities. After the collapse and virialization the system looked like a cluster of 
galaxies. Those early simulations of cluster formation, though producing cluster-like objects, 
signaled the first problem - simple model of initially isolated cloud (top-hat model) results 
in the density profile of the cluster which is way too steep (slope -3 -4) as compared with 
real clusters. The problem was addressed by Gunn & Gott [3], who introduced a notion of 
secondary infall in an effort to solve the problem. Another keystone work of those times is 
the paper by White [4], who studied collapse of 700 particles with different masses. It was 
shown that if one distributes the mass of a cluster to individual galaxies (=points), two-body 
relaxation will result in mass segregation which is not compatible with observed cluster. This 



was another manifestation of dark matter in clusters. This time it was shown that the dark 
matter of a cluster is not in galaxies. Survival of substructures in galaxy clusters was another 
problem addressed in the paper. It is quite remarkable that 20 years later we are still facing 
the same two problems. For numerical models (at least some of them) it is the excessive mass 
segregation, and, thus, a wrong structure of objects due to insufficient number of particles. 
For cluster dynamics it is the amount of substructure (buzzwords: "merging", "hydrostatic 
equilibrium" ) . 

First cosmological simulations (growth and collapse of initially small fluctuations in ex- 
panding Universe) were done at about the same time (mid 70s) . An example: Aarseth, Gott 
& Turner [5 ] studied the evolution of fluctuations using 5000 particles moving in expanding 
sphere. The particles were initially placed at random within the sphere (thus, the flat n = 
spectrum of initial fluctuations). The sphere was needed to handle boundary conditions: if a 
particle happens to hit the sphere it is reflected back like a ball in a billiard. At that time the 
correlation function completely dominated cosmological landscape. Simple hierarchical clus- 
tering arguments [6] indicated that in order to reproduce the observed slope of the correlation 
function of galaxies 7 = —1.77, the spectrum of initial fluctuations should be flat. Results of 
N-body simulations seems to confirm the conclusion. Another result, which followed from the 
same arguments was used to defeat a cosmological scenario proposed by Zeldovich [7 ] , which 
was called at that time "adiabatical fluctuations" and later resurfaced under the name of the 
Hot Dark Matter. Because the Zeldovich model suggested very distinct scale in the initial 
spectrum of fluctuations (almost no fluctuations on small scales), the argument was that 
it should not produce a scale-free power-law for the correlation function. As we know now, 
both arguments were wrong. And it was found thanks to numerical simulations. Efstathiou & 
Eastwood [8 ] using much better code (one of the first results with Particle-Particle-Particle- 
Mesh (P 3 M) code) with 20,000 particles showed that the model with the flat spectrum (initial 
random distribution) fails to produce power-law correlation function. Klypin & Shandarin [9] 
showed that the Hot-Dark-Matter (HDM) model predicts the power-law for the correlation 
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function. The HDM model is dead anyway (it does not have enough high-z objects), but the 
history of the model is really amazing: it was killed for sure so many times ... because of 
wrong reasons. 

Generation of initial condition with given amplitude and spectrum of fluctuations was a 
problem for some time. The only correctly simulated spectrum was the flat spectrum which 
was generated by randomly distributing particles. In order to generate fluctuations with 
power spectrum, say P(k) oc Aarseth, Gott & Turner [5] placed particles along rods. 
Formally, it generates the spectrum, but the distribution has nothing to do with cosmological 
fluctuations. As far as I know, the papers of Doroshkevich et al. [10] for two dimensions 
and Klypin & Shandarin [9] in three dimensions were the first where initial conditions were 
generated using the Zeldovich [7 ] approximation - the method which since then was used to 
generate initial conditions. 

Starting mid 80s the field of numerical simulations is blooming: new numerical techni- 
ques were invented, old ones were perfected, many publications (and, occasionally, results) 
are based on numerical modeling. To large extend, this have changed our way of doing 
cosmology. Instead of questionable assumptions and waving-hands arguments, we have tools 
of testing our hypothesis and models. As an example, I mention two analytical approximations 
which were validated by numerical simulations. The importance of both approximations is 
difficult to overestimate. The first is the Zeldovich approximation, which paved the way 
of understanding the large-scale structure of the galaxy distribution. The second is the 
Press-Schechter [11] approximation, which gives the number of objects formed at different 
scales at different epochs. Both approximations cannot be formally proved. The Zeldovich 
approximation formally is not applicable for hierarchical clustering. It must start with smooth 
perturbations (truncated spectrum). Nevertheless, numerical simulations have shown that 
even for the hierarchical clustering the approximation can be used with appropriate filtering 
of initial spectrum (e.g. [12,13,14]). The Press-Schechter approximation is also difficult to 
justify without numerical simulations. It operates with initial spectrum and linear theory, 
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but then (very long jump in assumptions) it predicts the number of objects at very nonlinear 
stage. Because it is not based on any realistic theory of nonlinear evolution (we just do not 
have one), it was an ingenious, but a wild guess. If anything, the approximation is based 
on a simple spherical top-hat model. But simulations show that objects do not form in 
this way - they are formed in a complicated fashion through multiple mergers and accretion 
along filaments. Still this very simple (and a very useful) prescription gives quite accurate 
predictions [15]. 

The following of this lecture is organized in the following way. Section 2 gives the equa- 
tions which we solve to follow the evolution of initially small fluctuations. A brief discussion 
of different methods is given in section 3. Effects of the resolution and some other technical 
details are also discussed in Section 3. Identification of halos ("galaxies") is discussed in 
Section 4. 

2. Equations of evolution of fluctuations in an expanding universe 

Usually the problem of the formation and dynamics of cosmological objects is formulated 
as N-body problem: for N point-like objects with given initial positions and velocities find 
their positions and velocities at any consequent moment. It should be remembered that this 
just a short-cut in our formulation - just to make things simple and avoid any discussion. 
While it still mathematically correct in many cases, it does not explain what we are doing. 
If we are literally to take this approach, we should follow the motion of zillions of axions, 
baryons, neutrinos, and whatever else our Universe is made of. So, what it has to do with the 
motion of those few millions of particles in our simulations? The correct approach is to start 
with the Vlasov equation coupled with the Poisson equation (and proper initial and boundary 
conditions). If neglect the baryonic component, which of course is very interesting, but would 
complicate our situation too much, the system is described by distribution functions /i(x, x, t) 
which should include all different clustered components i. For a simple CDM model we have 
only one component (axions or whatever it is). For more complicated Cold plus Hot Dark 
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Matter (CHDM) with few different types of neutrinos the system includes one DF for the 
cold component and one DF for each type of neutrino ([16]). The equations for the evolution 
of fi are: 

it + X fir " V4> lk = °' V2<P = 4nGa2 (P^ *) " = 4nGa 2 n dm 5 dm , 

5 dm (x,t) = (pdm- (Pdm))/(pdm)), /0 dm (x, t) = a" 3 ^ m i / d 3 x/i(x, A, t). 

Here a is the expansion parameter, Odm is the contribution of the clustered dark matter to 
the mean density of the Universe, mi is the mass of a particle of ith component of the dark 
matter. The solution of the Vlasov equation can be given in terms of characteristic equations, 
which look like equations of particle motion. The distribution function is constant along each 
characteristic. Complete set of characteristics is equivalent to the Vlasov equation. Of course, 
we can not have the complete set, but we can follow the evolution of the system (with some 
accuracy) if we select a representative sample of characteristics. One way of doing this would 
be to split initial phase space into small domains, take only one characteristic as representative 
for this volume element, and follow the evolution of the system of the "particles" in self- 
consistent way. In models with one "cold" component of clustering dark matter (like CDM 
or ACDM) the initial velocity is a unique function of coordinates (only "Zeldovich" part is 
present, no thermal velocities). This means that we need to split only coordinate space, not 
velocity space. For complicated models with significant thermal component (CHDM), the 
distribution in full phase space should be taken into account. Usually, this is simulated by 
placing 1-10 particles in each coordinate with velocities, which mimic the distribution of real 
"thermal" velocities. Depending on what we are interested in, we might split initial space 
into equal-size boxes (typical setup for PM or P 3 M simulations) or we could divide some area 
of interest (say, where a cluster will form) into smaller boxes, and use much bigger boxes 
outside the area (to mimic gravitational forces of the outside material) . In any case, the mass 
assigned to a "particle" is equal to the mass of the domain it represents. Now we can think 
of the "particle" either as a small box, which moves with the flow, but does not change its 



original shape, or as a point-like particle. Both presentations are used in simulations. None 
is superior to another. 

There are different forms of final equations. If we choose "momentum" p = a 2 x as 
effective velocity (v pec = p/a,) and change the independent variable from time t to expansion 
parameter a, then the equations are 

dp V0 rfx p 2 2 

3- = — , 3- = -r-y, V (p = AnGiloa d dm , 

^ da a da aa z 
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® ~ y^O + ^curv,O a + ^A,O a3 7 ^0 + ^curv,0 + ^A,0 — 1 

where Qq, ^curv,0) an d ^a,o are the density of the matter, effective densities of the curvature 
and cosmological constant in units of the critical density at z = 0. The curvature contribution 
is positive for negative curvature. 

3. Methods 

There are many different numerical techniques to follow the evolution of a system of many 
particles. For earlier reviews see Hockney & Eastwood [17] and Sellwood [18]. The most 
frequently used methods for cosmological applications fall in three classes: Particle Mesh 
(PM) codes, Particle-Particle/Particle- Mesh (P 3 M) codes, and TREE codes. All methods 
have their advantages and disadvantages. 

PM code. It uses a mesh to produce density and potential. As the result, its resolution 
is limited by the size of the mesh. Largest simulations were done by the author: 800 3 mesh 
with 3 x 256 3 = 1.5 x 10 8 particles. New parallel supercomputer SP2 at Cornell will be used 
to run simulations with 1600 3 = 4.096 x 10 9 mesh. There are two advantages of the method: 
i) it is fast (the smallest number of operations per particle per time step of all the other 
methods), ii) it typically uses very large number of particles. The later can be crucial for 
some applications. There are few modifications of the code. There are few variants of PM 
code, "plain- vanilla" PM was described by Hockney & Eastwood [17]. It includes Cloud-In- 
Cell density assignment and 7-point discrete analog of the laplacian operator. Higher order 
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approximations improve the accuracy on large distances, but degrade the resolution (e.g. [19] 
). In an effort to reduce the order of approximation and to increase the resolution, Melott [20] 
introduced staggered mesh. It gives a better resolution on cell-size distances, but particles get 
self-forces (an isolated particle experiences a force from itself), which might be not a welcome 
feature. No serious testing of the method was done. 

P 3 M code is described in detail in [17] and [21 ]. It has two parts: PM part, which takes 
care of large-scale forces, and PP part, which adds small-scale particle-particle contribution. 
The simulations usually have 64 3 -100 3 particles. The best known to me simulations were 
done by Ma & Bertschinger [22] for CHDM model: 128 3 of cold particles and ten times 
more hot particles. Because of strong clustering at late stages of evolution, PP part becomes 
prohibitively expensive once large objects start to form in large numbers. Significant speed 
is achieved in modified version of the code, which introduces subgrids (next levels of PM) in 
areas with high density [23]. With modification the code runs as fas as TREE code even for 
heavily clustered configurations [23]. 

TREE code is the most flexible code in the sense of the choice of boundary conditions 
[24 , 25 , 26 . It is also more expensive than PM: it takes 10-50 times more operations. 
Bouchet & Hernquist [27] and Hernquist, Suto & Bouchet [28] extended the code for the 
periodical boundary conditions, which is important for simulating large-scale fluctuations. 
An interesting blend of PM and TREE code as done by Xu [29]. The code might be easier 
than usual to implement on parallel supercomputers. 

Multigrid methods were introduced long ago, but only recently they started to show a 
potential to produce real results [30, 31 , 32]. It worth of paying attention if a "multigrid" 
code is really a fully adaptive multigrid code. Some of the them are actually two-level mesh 
codes. Those codes provide increased resolution (by factor 4-8) inside one predefined region. 

What code is the best? Which one to choose? Should we through away results of 
low resolution codes? (Tempting, but I would not do that). There is no unique answer - 
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everything depends on the problem, which we are addressing. For example, if we are interested 
in explanation of the large-scale structure (filaments, voids, Zeldovich approximation, and so 
on), PM code with 256 3 mesh is sufficient. It takes only one night to make a simulation on a 
(good) workstation. There is a very long list of problems like that. But if you intent to look 
for galaxies in the large-scale environment, much better resolution is needed. Sometimes a 
minimum dynamical range of 10 4 is quoted: we need to go to 10 kpc to see a galaxy, and 
we need to go up to at least 100 Mpc in order to have large voids and superclusters. What 
we should also include in this "wish-list" is a required mass range. If we would like to have 
something, which looks like a 10 11 M & galaxy, it should consist of many particles. Say, 
30-100 particles. This gives 10 9 M© per particle and 410 3 particles for 100 Mpc box. All 
existing methods fall short of these requirements: PM is lacking resolution (a factor of ten), 
others do not have enough mass resolution (a factor of ten). 

4. Distribution of the matter and effects of resolution 

Figure 1 shows an example of structures observed in cosmological models. The following 
model with cosmological constant was chosen: fi = 0.30, h = 0.70, as = 1.1, <5rms-ps = 
2l.8fiK. The simulation was done with the PM code with 256 3 particles, 800 3 mesh, Aa = 
0.003, box size 50/i _1 Mpc. Dark matter particles in a l/i _1 Mpc slice are shown. Different 
structures can be identified on the plot. Large empty (of galaxies) voids, long filaments, few 
groups of galaxies, and many isolated galaxies. While it is not clear from this plot whether 
we actually see filaments or sheets, Kofman et al. [33] recently argued that they actually 
are filaments. It was a long-standing problem: the Zeldovich approximation was predicting 
"pancakes" and what was "numerically observed" looked more like filaments. The argument 
from the Zeldovich approximation is still valid. If we smooth initial density field and then 
randomly pick up a point at very early stage and trace its trajectory, then the particle will 
first hit a sheet, not a filament, because its middle eigenvalue of the deformation tensor is 
zero (on average) and the smallest one is negative. (Positive sign of an eigenvalue implies 
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collapse along corresponding direction. Only one positive eigenvalue means collapse along 
only one directions - a sheet forms). But Kofman et al. argue that this is not we are doing. 
We observe particle at a given moment. In this case we are asking what is the shape of 
an isodensity surface with density of few mean densities. This condition results in positive 
middle eigenvalue while keeping the expansion along the third axis - this is a filament. 

In order to discuss objects on smaller scales we need to understand what happens with an 
isolated object, which size is close to the limit of resolution. In all codes the resolution (either 
a cell size or a softening parameter) is constant in comoving not proper coordinates. This 
means that the resolution actually decreases with time. (Look at the situation from a positive 
angle - the resolution gets much better if we go to high redshifts). For example, our resolution 
is 50 kpc and we have a galaxy with radius 50 kpc collapsing at z = 4, a = 0.2. The resolution 
at that redshift was 10 kpc - enough to reasonably resolve the object. Now, our resolution 
drops very severely. If this happened quickly (on dynamical scale), the galaxy would explode: 
its total energy would be positive. Fortunately, this does not happen because the time-scale 
of growth of the softening parameter is equal to the expansion rate of the Universe, and the 
later is always longer than the dynamical time. Thus the galaxy will adiabatically expand 
once the softening length gets close to its radius. Figure 2 illustrates this behavior. A system 
of 300 particles with initial velocities slightly below virial velocities was simulated using a 
PP code. Initially the particles were randomly distributed in a sphere with 100 kpc radius. 
Parameters were scaled in such a way that the total mass of the system is 10 12 M , initial 
moment corresponds to a = 0.2, z = 4. Collapse happens at a = 0.25, z = 3. The radius 
of half mass is 50 kpc after collapse. The system was run once with softening parameter 
e = 10 kpc in proper coordinates, and another time with e = 10 kpc in comoving coordinates. 
Figure 2 shows snapshots of the simulations. Comoving coordinates are chosen to display the 
particles. Coordinate scale is in units of 100 kpc. Bottom row shows the evolution of the 
system simulated with constant proper resolution, while the top row is for constant comoving 
resolution. The circles show the size of resolution. Radius of each circle is 1.5e. At this radius 
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the error in force at 1 cell radius in PM code is equal to the error in PP code. As time goes 
on the system shrinks in comoving coordinates. But once its size is close to the resolution, 
the shrinking stops. 

This just confirms what is found in numerical simulations: isolated unresolved "galaxies" 
are small, almost spherical balls with 2 grid cells across their diameters. 

In order get a better insight on the effects of the resolution and to understand what 
should be expected if we significantly increase the resolution, two simulations were run with 
exactly the same initial conditions, but with twice different resolution. The same cosmological 
model as before was chosen (ACDM with fi = 0.3, h = 0.70). PM code with 128 3 particles, 
mesh 450 3 and 216 3 , box size 20/i _1 Mpc was used. This gives the resolution of 44/i _1 kpc 
and 92/i _1 kpc for the two runs and mass for a particle mi = 3.1 x 10 8 /i _1 M . Figure 
3 shows particles in a small 3x3 Mpc window in the simulations at z = 3. Only particles 
with estimated density indicated in the figure are shown and only 1/4 of all particles is 
displayed. Because the same isolated clump would have smaller radius and, thus, higher 
density with twice better resolution, density limits were slightly adjusted to take into account 
the difference in the resolution. Figure 4 shows the same window at z = 0. While plots 
for higher resolution show more small clumps, all large clumps are found in both plots. The 
differences are as expected: higher resolution results in more compact and dense objects. The 
most significant result lies in what was feared, but not found. The largest object in z = 
plot has radius of about 300 kpc. The usual question was: having a low-resolution run with 
100 kpc resolution, how do you know that with better resolution it will not split into, say 
into two objects? Our results indicate that it does not happen: with more then twice better 
resolution the object got a bit smaller, a bit rounder, got few tiny satellites, which were barely 
seeing with the low resolution, but it does not show any tendency to break int large peaces. 
Figures 5 and 6 show another examples of comparison of the two simulations. In this case 
we deal with a group of galaxies. The result is just the same: more small objects, the same 
large objects, no tendency for splitting of large halos into smaller ones. 
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While visual analysis indicates that we might be quite ok with relatively low resolution 
if only large halos (above 10 10 M & ) are considered, it is not that easy to make quantitative 
analysis. Comparison of coordinates, velocities, and densities of individual particles in both 
simulations would be a natural test, but it just fails. The reason is the divergence of trajec- 
tories in dynamical systems: small differences in initial coordinates result in large differences 
after few dynamical times. The winding problem of the spiral pattern is just one of the man- 
ifestation of this divergence. As the result, deviations of coordinates X44 — x 92 and velocities 
Kr,44 — Vx,92 as functions of density (Figures 7 and 8) indicate mainly sizes of virialized object, 
which exist at different redshift. Even in spite of the divergence of the trajectories, extremely 
large fraction of particles - 99% - indicate the difference in coordinates less than 2 cell sizes 
and difference in velocities less than 100 km/s. 

5. Halo identification and overmerging problem 

There are different methods of identifying collapsed objects (halos) in numerical simula- 
tions. 

Friends-Of-Friends (FOF) algorithm was used a lot and still has its adepts. If we 
imagine that each particle is surrounded by a sphere of radius bd/2, then every connected 
group of particles is identified as a halo. Here d is the mean distance between particles, and 
b is called linking parameter, which typically is 0.2. Dependance of groups on b is extremely 
strong. The method stems from an old idea to use percolation theory to discriminate between 
cosmological models. Because of that, FOF is also called percolation method, which is wrong 
because the percolation is about groups spanning the whole box, not collapsed and compact 
objects. FOF was criticized for failing to find separate groups in cases when those groups 
were obviously present [19]. The problem originates from the tendency of FOF to "percolate" 
through bridges connecting interacting galaxies or galaxies in high density backgrounds. 

DENMAX tried to overcome the problems of FOF by dealing with density maxima 
[19 ,34]. It finds maxima of density and then tries to identify particles, which belong to each 
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maximum (halo). The procedure is quite complicated. First, density field is constructed. 
Second, the density (with negative sign) is treated as potential in which particles start to 
move as in a viscous fluid. Eventially, particles sink at bottoms of the potential (which are 
also maxima density). Third, only particles with negative energy (relative to their group) 
are retained. Just as in the case of FOF, we can easily imagine situations when (this time) 
DENMAX should fail. For example, two colliding galaxies in a cluster of galaxies. Because of 
large relative velocity they should just pass each other. In the moment of collision DENMAX 
ceases to "see" both galaxies because all particle have positive energies. That is probably 
quite unlikely situation. The method is definitely one of the best at present. The only problem 
is that it seems to be too complicated for present state of simulations. 

"Overdensity 200". There is no name for the method, but it is often used. Find 
density maximum, place a sphere and find radius, within which the sphere has the mean 
overdensity 200 (or 177 if you really want to follow the top-hat model of nonlinear collapse). 
In Figure 5 all particles in the central ~1 Mpc area will be just one halo. The mass function 
of dark halos, constructed in this way, has a very long tail extending into masses, which are 
far too large for individual galaxies. The existence of this tail is often called "overmerging 
problem". Figure 9 gives comparison of mass functions of dark halos n(> M) in previous 
two simulations with different resolutions. The full curve shows the mass function for high 
resolution run. Note that high mass tail of the mass functions does not depend on resolution, 
which indicates convergence of the results. The mass function flattens in both simulations at 
different level (more small-mass objects with high resolution), but at the same mass, which 
corresponds to 20-30 particles per object. 

Actually, there are two "overmerging" problems. One arises because of naive application 
of the top-hat model. Consider our Galaxy and Andromeda Nebula as a testbed. If we assume 
that both galaxies have flat rotation curves at least up to the radius of 100 kpc (which might 
be true), then overdensity at 100 kpc is about 4000 for our Galaxy and about twice that for 
Andromeda Nebula. The overdensity at half way between galaxies is about 200-300. Thus, 
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if we had fantastic simulation with our Galaxy and M31 as they are in the Universe, our 
"overdensity 200" method would lump them together. There are different ways of fixing the 
problem. DENMAX probably would find both galaxies. 

Unfortunately, there is real "overmerger" problem, which cannot be avoided that easily. 
For a long time it was assumed that appearance of extremely large halos ( "overmergers" ) is 
due to insufficient resolution. Numerical simulations indicated that with smaller and smaller 
mass of each particle, and with better and better resolution we see more and more substructure 
in overmergers. But it seems that the tendency have stopped. Now extra resolution does not 
split large overmergers. The problem is real and it is not a numerical artifact, which can 
be cured by much better resolution. Figures 3-6 illustrate this result. Recently, Moore et 
al. [35] and van Kampen [36] gave a very plausible explanation: tidal forces. If a small 
halo ("galaxy") falls into a larger one ("group"), it will be tidally disrupted at a distance 
where its central density is equal to the mean density of the large halo at that distance (for 
"group" with p oc r -2 ). Because the core of the large halo likely was formed earlier than the 
small halo, it has higher density. The only chance for the small halo to survive is to stay in 
peripheral parts of the group. Numerical simulations show that this is what happens. 
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Figure Captions 

Figure 1 An example of structures observed in cosmological models with cosmological con- 
stant: fio = 0.30, h = 0.70, ag = 1.1, Qrms-ps = 21. 8 fj,K. The simulation was done with 
the PM code with 256 3 particles, 800 3 mesh. Dark matter particles in a l/i _1 Mpc slice are 
shown. 

Figure 2 Evolution of a system of 300 particles in comoving coordinates. Initially the 
particles were randomly distributed in a sphere with 100 kpc radius. Collapse happens at 
a = 0.25, z = 3. The system was run once with softening parameter e = 10 kpc in proper 
coordinates (bottom row), and another time with e = 10 kpc in comoving coordinates (top 
row). Coordinate scale is in units of 100 kpc. The circles show the size of resolution. Radius 
of each circle is 1.5e. 
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Figure 3 Particles in a small 3x3 Mpc window at z = 3 in ACDM model with Q = 0.3, h = 
0.70. Only particles with estimated density indicated in the figure are shown and only 1/4 of 
all particles is displayed. 
Figure 4 The same window at z = 0. 

Figure 5 An example of a group of galaxies simulated with with different resolutions. 

Figure 6 The sane as in Figure 5, but with different density threshold. 

Figure 7 Deviations of coordinates and velocities in high and low resolution simulations at 

z = 3. 

Figure 8 The same as in Figure 7, but at z = 0. 

Figure 9 Comparison of mass functions of dark halos n(> M) in two simulations with 
different reolutions (44/i _1 kpc and 92/i _1 kpc). The full curve shows the mass function for 
high resolution run. Note that high mass tail of the mass functions does not depend on 
resolution, which indicates convergence of the results. 
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